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BARRIER METHODS FOR CRITICAL EXPONENT PROBLEMS IN 
GEOMETRIC ANALYSIS AND MATHEMATICAL PHYSICS 

JENNIFER ERWAY AND MICHAEL HOLST 

Abstract. We consider the design and analysis of numerical methods for approxi- 
mating positive solutions to nonlinear geometric elliptic partial differential equations 
containing critical exponents. This class of problems includes the Yamabe problem 
and the Einstein constraint equations, which simultaneously contain several challeng- 
ing features: high spatial dimension rt ^ 3, varying (potentially non-smooth) coef- 
ficients, critical (even super-critical) nonlinearity, non-monotone nonlinearity (arising 
from a non-convex energy), and spatial domains that are typically Riemannian mani- 
^1 folds rather than simply open sets in _R". These problems may exhibit multiple solu- 

I tions, although only positive solutions typically have meaning. This creates additional 

^ complexities in both the theory and numerical treatment of such problems, as this fea- 

^^ ture introduces both non-uniqueness as well as the need to incorporate an inequality 

Cn constraint into the formulation. In this work, we consider numerical methods based on 

^^ Galerkin-type discretization, covering any standard bases construction (finite element, 

(~| spectral, or wavelet), and the combination of a barrier method for nonconvex optimiza- 

Qh tion and global inexact Newton-type methods for dealing with nonconvexity and the 

(— I presence of inequality constraints. We first give an overview of barrier methods in non- 

+-^ convex optimization, and then develop and analyze both a primal barrier energy method 

a for this class of problems. We then consider a sequence of numerical experiments using 

this type of barrier method, based on a particular Galerkin method, namely the piecewise 
linear finite element method, leverage the FETK modeling package. We illustrate the be- 
T-H havior of the primal barrier energy method for several examples, including the Yamabe 

K^ problem and the Hamiltoruan constraint. 
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1. Introduction 

In this article we consider the design and analysis of numerical methods for approx- 
imating positive solutions to nonlinear geometric elliptic partial differential equations 
containing critical exponents. These types of problems arise regularly in geometric anal- 
ysis and mathematical physics, examples of which include the Yamabe problem and 
the Einstein constraint equations [8, 9]. These problems often simultaneously contain 
several challenging features, including spatial dimension n ^ 3, varying and potentially 
non-smooth coefficients, critical (or even super-critical) nonlinearity, non-monotone non- 
linearity (arising from a nonconvex energy), and spatial domains that are typically Rie- 
mannian manifolds rather than simply open sets in i?". For these types of problems, there 
may be multiple solutions, although only positive solutions typically have mathematical 
and physical meaning. This creates additional complexities in both the theory and numer- 
ical treatment of such problems, as this feature introduces both non-uniqueness as well 
as the need to incorporate an inequality constraint into the formulation. In this work, we 
consider numerical methods based on Galerkin-type discretization, covering any stan- 
dard bases construction (finite element, spectral, or wavelet), and the combination of a 
barrier method for nonconvex optimization and global inexact Newton-type methods for 
dealing with nonconvexity and the presence of inequality constraints. Our goal is to de- 
velop reliable methods for computing positive approximate solutions to these types of 
nonlinear problems. 

Critical exponent problems arise in a fundamental way throughout geometric analysis 
and mathematical general relativity. One of the seminal problems in this area is the 
Yamabe Problem: Find u E X such that 

~8AgU + Ru = RuU^ in 11, (1.1) 

u>0, (1.2) 

where 17 is a Riemannian 3-manifold, g is the positive definite metric on f2, A^, is the 
Laplace-Beltrami operator generated by g, R is the scalar curvature of g, and Ru is the 
scalar curvature corresponding to the conformally transformed metric: 

g = (f)^g. (1.3) 

The coefficients R and Ru can take any sign. The Banach space X containing the solu- 
tion is an appropriate Sobolev class ^^'^'^(17) for suitably chosen exponents s and p. If 
the manifold 17 has a boundary, then boundary conditions are also prescribed, such as 
M = 1 on an exterior boundary to 17. In the case that 17 C M^, and gij = 5ij, then A^, re- 
duces to just the Laplace operator on 17. With the presence of the term u^ and the spatial 
dimension being three, this is an example of a critical exponent problem; such problems 
are known to be difficult to analyze as well as to simulate numerically. The presence of 
the inequality constraint (only positive solutions have mathematical and physical mean- 
ing) creates additional complexities in both the theory and numerical treatment of such 
problems. Prior work on numerical methods for critical exponent semilinear problems 
has focused primarily on the development of adaptive methods for recovering solution 
blowup; cf. [3, 2]. 

Outline of the paper The structure of the remainder of the paper is as follows. In §2, 
we give a more detailed overview of the class of geometric PDF problems of interest, 
including both the Yamabe problem and the Hamiltonian constraint in the Finstein equa- 
tions. As part of the discussion, we derive the linearized Hamiltonian constraint, and 
construct an artificial nonconvex "energy" functional, which gives rise to the Hamilton- 
ian constraint as a condition for its stationarity. In §3, we give an overview of barrier 
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methods in nonconvex optimization. In §4 we then develop and analyze a primal bar- 
rier energy method for this class of problems. Finally, in §5 we consider a sequence of 
numerical experiments using this type of barrier method, based on a particular Galerkin 
method, namely the piecewise linear finite element method, leverage the FETK model- 
ing package. We illustrate the behavior of the primal barrier energy method for several 
examples, including the Yamabe problem and the Hamiltonian constraint. We draw some 
conclusions in §6. 

2. Elliptic Problems in Geometric Analysis and Relativity 

While one of our motivations here is to develop methods for the Yamabe problem and 
similar problems arising in geometric analysis, we are also interested in a related, more 
general problem arising in general relativity. The Einstein equations, which represents 
Einstein's 1915 theory of gravity, are a coupled hyperbolic-elliptic system that governs 
the deformation of the underlying metric of spacetime in response to the distribution and 
dynamics of matter and energy density. The elliptic part of the system, known as the 
Einstein constraint equations, or the coupled Hamiltonian and momentum constraints, 
are of great interest in both the mathematical and numerical relativity research commu- 
nities. This elliptic system must be satisfied by initial data used to evolve the metric 
forward in time with the hyperbolic portion of the Einstein equations (called the evolu- 
tion equations), and the constraints must also be satisfied at all points in time during the 
evolution. 

The Einstein constraints have all of the difficult features of the Yamabe problem, plus 
more: three spatial dimensions, non-flat Riemannian manifold spatial domain, critical ex- 
ponent, non-monotone nonlinearity, negative exponent powers (non-polynomial rational 
nonlinearity, giving rise to singularities at the origin), possibly non-smooth coefficients, 
possible non-uniqueness, physical positivity requirement, and the structure of an elliptic 
system for two variables, with lack of a variational structure (the equations do not arise 
as the Euler condition for stationarity of an underlying energy functional). However, one 
of the difficulties is not present in an important physical situation known as the constant 
mean curvature (CMC) case: in this situation, the coupled elliptic system of the Hamil- 
tonian and momentum constraints decouple into the separate constraints, both of which 
have separate variational structure. However, all of the other difficulties remain, and the 
Hamiltonian constraint alone may be viewed as a generalization of the Yamabe problem. 
An overview of the Einstein constraints, including the CMC case, can be found in [8, 9]. 
Here we will consider only the CMC case, and focus on the Hamiltonian constraint, also 
known as the Lichnerovich equation. 

p 2 2 

-V ■{aVu) + —u = -—u^ + —u''^ + 2'Kpu'^ in fi, (2.1) 

8 12 8 

{aVu{x,y,z)) ■ n + cu = qn on 9^^^, (2.2) 

u = gn on S/jfi, (2.3) 

where n is the unit normal and dVl = d^^ [j d^^ and do^ fl ^a^^ — ^- Here R{x), 
t'^{x), (T^(a;) : f] C 5R^ -)■ 3fi and p, a^, r^ > for aU x e il. Also it assumed that 
there exist positive constants Ci and C2 such that — Ci < R{x) < C2. Also, r may 
be considered to be constant so that CMC decoupling occurs. Reasonable values for a 
are such that > cr^ < C3 where C3 may be as large as 10^. The strong form of the 
constraints given in (2.1-2.3) can be transformed into the weak form, reformulating the 
problem using fewer derivatives. 
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2.1. Weak Formulation. The weak formulation is obtained by taking the L^-inner prod- 
uct over fi with all test functions i; G -f^o(fi) = {u e H^{Q) : t; = on Q^} and (2.1), 
yielding: 



-aAu 



R 



-\. ^\.-7 



u H u u — 2'Kpu ] V dx = 0. 

12 8 



Green's first identity states: 

/ (V ■ z)v dx = {n ■ z)v ds — I z ■ Vf dx. 
Jn Jan Jn 

Taking z := aVu and recalling that t; = on doi^), we obtain: 

/ (aVu) ■ Vv dx — [n ■ dVu)v ds + k{u)v dx = 0, 



(2.4) 



where k{u) 



g-M + -^2^ g It 



2tx pu-^. Using (2.2) in (2.4), yields: 



aVn) ■ Vf + k{u)v dx — I {g^ — cu) v ds = 0. 
n JdM^ 



Thus, the weak form of (2.1-2.3) is given by: 

Findw G X = H}){n) n [u.,u+] s.t. {f{u),v) = 0, for all w G H^iCl), 
where //^(fi) = {u e H^{Q) : u = g^ on ^d} and 

{f{u),v)= / [{aVu)-Vv + k{u)v] dx + {cu-gN)v. 



(2.5) 



(2.6) 



Note that we have constructed the space X = Hl){Q,) fl [«_,«+] in which to look for 
solutions to the problem based on the need to "guard" the nonlinearity from blow-up at 
the origin. This potential blowup is due to the negative powers appearing as part of the 
non-polynomial, rational form of the nonlinearity. The pointwise interval [n_ , u+] , which 
can be strictly negative or strictly positive, can be shown to contain one or more solutions 
using maximum principles and fixed-point arguments; cf. [8, 9]. The numerical methods 
we develop later in the paper will incorporate this type of "guarding" in the discrete 
formulation. 

The Gateaux derivative (Df{u)w, v) of the weak nonlinear form (2.6) is needed for 
use in Newton-like algorithms. It is computed formally as 



di 



f{u + tw) {v) 



t=0 



d_ 
di 



+ 



d_ 
Ft 



/ [{aV{u + tw)) -Vv + k{u + tw)v] dx 
Jq 

(c (m + tw) — g^) V ds 



t=o 



divfl 



t=0 



(2.7) 



giving that 



{Df{u)w,v) 



/ [{aVw) ■ Vf + k'{u)wv] dx + cwv ds, (2.8) 

Jq JdNO, 



where k'{u)wv = ^wv + jkt'^u^wv + Ict^m ^wv + Qnpu ^wv. 
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2.2. The energy. The weak formulation can be viewed as a zero-finding problem; al- 
ternatively, it may be viewed as the problem of finding a critical point of an energy 
functional. 

Theorem 2.1. u is a solution to the weak form (2.5) if and only ifu is a critical point of 
the energy functional 



J{u) 



'aVu)-Vu+-Ru^ ' ^-2.,6, 1 „2.,-6 , ^..,-2 



1 

16' 



+ ^^ w" + ^^^" + ^P« 



+ - 



cu ds 



9jvO 



dx 
qnu ds. (2.9) 



(9jvf! 



Proof We prove this theorem by showing that Gateaux derivative of J{u) is exactly the 
weak formulation (2.5). We have 



d_ 
di 



J{u + tv) 



t=0 



d_ 
di 



Q 



- (aV(u + tv)) ■ V(m + tv) + —R{u + tv)^ 



+— r^(w + tv)'^ + —(y^{u + tv)-^ + Txp{u + tv)-"^ 
1 2, 4o 



dx 



t=o 



+ 



di 



9jvn 



1 2 

-c{u + tv) ~gN{u + tv) 



ds 



t=o 



p 2 2 

(aVw) ■ Vf H uv H u^v u^'^v — 2tx pu^^v 
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+ / {cu — gN)v ds 
U\u),v). 



dx 



U 



Theorem 2.1 can be restated in terms of stationarity. 

Corollary 2.2. The energy functional J{u) in (2.9) is stationary at u if and only ifu is a 
solution to the weak formulation (2.5). 

From an optimization standpoint, it is reasonable to ask under what conditions can 
J(m) be minimized to obtain the solution to the weak formulation. Certainly, if J(m) is a 
convex energy functional then a weak solution can be found by minimizing J{u). Con- 
vexity implies uniqueness of solutions, which in this application is often not the case. 
Nevertheless, Newton's method may be used to find a critical point of the energy, requir- 
ing the need for a second derivative. From Theorem 2.1, the second Gateaux derivative 
of the energy functional J{u) is exactly D{f{u),v) (see (2.8)). 

Example 1. Set a = 1.0, R = 1.0, r = 0.1, a = 0.2, p = 0.1, c = 1.0, and 
gj\f = —1.0 with Vl chosen to be a single hole domain at the origin in three-dimensional 
space with only Robin boundary conditions. In this simple case, since R > and c > 0, 
the energy functional J(m) is convex on the domain of positive functions u > 0, and 
thus, minimizing the energy functional J over n > is equivalent to solving the weak 
formulation. Similarly, if i? > and c > 0, then J{u) is also convex over m < 0- 
allowing for the existence of both strictly negative and positive solutions. This yields a 
convex energy functional J with one positive solution m > and one negative solution 
u < (see Section 5). 
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Figure 1 . A plot of the integrand /(m) with u = x, a = 2, R 

T = V72, a = 748, p = I/tt, and fi = [0.4, 3]. 



-1000, 



Example 2. Set a = 2.0, r = v72, a = v48, and p = l/vr, and let Q be the 
one-dimensional closed subset [0.1, 10]. With these choices, the energy reduces to: 



J{u) 



10 



0.1 



(Vm) ■ Vu + 



16 



Ru"^ + u^ + u~ 



+ u 



dx. 



(2.10) 



Let I{u) denote the integrand in (2.10). For u = x, the second derivative of I{u) is given 
by 

dx'^ 8 
Thus, if i? < and sufficiently negative, / has an inflection point; otherwise, / is a 
strictly positive, convex function of u. For any value of R, I{u) is an even function since 
I{u) = I{—u). The integrand I{u) with R = —1000 is plotted in Figure 1 in coordinate 
pairs (x, /(«)) with u = x on the interval [0.4, 3]. The second derivative of I{u) confirms 
this is a nonconvex function of u. 

When the energy functional is convex, a weak solution can be found by minimizing 
the energy functional using a Newton-like iteration. However, when the energy is non- 
convex, a weak solution may be a maximizer or saddlepoint of the energy. In this case, 
Newton's method can be used to find a stationary point of the energy functional with 
progress towards a stationary point being guaranteed by computing steps along the New- 
ton direction that yield sufficient decrease in a merit function. The positivity constraint 
u > can naturally be enforced when solving the Lichnerovich equation (2.1-2.3) using 
a safe-guarded Newton method. Unlike the Lichnerovich equation, the Yamabe Problem 
(1.1) does not have a singularity at u = 0; thus, a barrier method can be used to help 
enforce the positive inequality constraint on u. In the following section, we consider 
barrier methods for nonconvex optimization and develop a primal barrier energy method 
for this class of problems. 



3. Barrier methods for nonconvex optimization 

Barrier methods are the most widely-used type of interior method for general nonlinear 
inequality-constrained optimization problems of the form: 



subject to c{x) > 0, 



mminiize 



(3.1) 
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where / : 3f?" — t- 5R is assumed to be twice-continuously differentiable and c : 5R" — > SfJ™ 
is an m-vector of constraints. Generally speaking, barrier methods seek to minimize a 
composite function that both resembles the original function and naturally prevents in- 
feasible iterates. Today, the most widely used barrier function is the classical logarithmic 
barrier function: 

m 

B^{x) = f{x)- ^i'^\nci{x). (3.2) 

i=l 

Notice that if /i is small, the barrier function resembles the original function; moreover, 
this function inherits the smoothness associated with the original problem but is only 
defined in the strict interior of the feasible region for the original problem (3.1). 

The classical barrier method solves (3.1) by minimizing B^{x) for a decreasing se- 
quence of positive jj,. Given a fixed /x > 0, first-order optimality conditions for x* to be 
a minimizer of the barrier function is that Vi?^(x*) = 0, i.e., 

V/(x*) - J{x*Yy{x*) = 0, (3.3) 

where J denotes the constraint Jacobian and yi{x) = n/ci{x). Alternatively, y is inter- 
preted as a vector of Lagrange multipliers associated with the original inequality problem 
(3.1). Moreover, Ci{x)yi = jj, can be viewed a perturbation of the complementarity con- 
dition for a first-order KKT point. 

When c(x) = x, the Newton equations for minimizing B^{x) are given by 

(VV(x) + /idiag(X-2)) p=- [V/(x) - /iX-^] , (3.4) 

where X~^ = [l/x{, I/X2, • . . , l/a;{]^ for j G {1, 2}, and diag(x) is the diagonal matrix 
whose 2th diagonal entry is Xj. To ensure global convergence, a line search must be used 
to satisfy a sufficient decrease criteria (e.g., the Armijo or Wolfe conditions). Once a 
minimizer x(/i) of B^{x) is computed, jj, is reduced and the process is repeated. Notice 
that the subsequent minimization can be warm-started by choosing the minimizer of the 
previous barrier function as the initial guess. Algorithm 3.1 summarizes the classical bar- 
rier method for solving (3.1) with c(x) = x and using the Amijo condition for sufficient 
descent. 



Algorithm 3.1: Classical Barrier Method. 

Choose a;o > 0, yU > 0, r^ G (O, ^), and 7 G (0,1); 
Set /c = 0; 

while not converged do 

Compute x(/i), an unconstrained minimizer of B^{x): 
while not converged do 
Solve (3.4) to obtain p; 
Use the "99% rule" to compute amax; 

Compute a G (0, ctmax] such that B^{x + ap) < B^{x) + -qaV B ^{xY p; 
X ■k- X + ap; 
end do 
Xk+i ^ x{fi); 

k ^k + 1; 
end do 
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When computing a step length, it is necessary to safeguard the step to avoid taking 
a step into the infeasible region. For linear constraints, the so-called "99%" rule may 
be invoked (e.g., see [6]) that states that if Axi < for at least one i: (i) Compute 
dmax = max{— (xj/Axj) : Axj < 0}; (ii) Set the maximum step length to be a = 
max{0.99 x amax, !}• 

We now state two theorems that summarize the convergence of the classical barrier 
method. The first theorem governs the local convergence of a sequence of minimizers 
of the classical barrier function. Before stating this theorem, we state the following 
theorem, which will be used in the proof of the theorem on local convergence. 

Theorem 3.1. Consider problem (3.1), where / : 3?" — )• 3f? and c : 3fi" — )■ W^ are 
continuous. Let M denote the set of all local constrained minimizers with objective 
function value f*, and assume that f* has been chosen so that M is nonempty. Assume 
further that the set M* ^ Af is a nonempty compact isolated subset of M. Then there 
exists a compact set S such that M* lies in int(5) fl J-" and for any feasible point x in 
S but not in M*, f{x) > /*. Furthermore, every point x* in A/"* has the property that 
f{x*) = f* = mm f (x) for all X G 5 fl J". 

Proof. See [5, Theorem 7] and [12, Theorem 6]. D 

The following theorem and proof regarding convergence of the minimizers of the bar- 
rier function is found in [6]. 

Theorem 3.2. Consider problem (3.1), where / : 3fJ" — > 3f? and c : 3f?" — )■ SfJ™ are 
continuous. Let T denote the feasible region, let M denote the set of minimizers with 
objective function value /*, and assume that M is nonempty. Let {/ifc} be a strictly 
decreasing sequence of positive barrier parameters such that limfc__>oo l^k = 0. Assume 
that 

(a) there exists a nonempty compact set M* of local minimizers that is an isolated 
subset ofAf; 

(b) at least one point in Af* is in the closure o/ strict (J-"). 

Then the following results hold: 

(i) there exists a compact set S such that A/"* C int (S) fl J-" and such that, for any 

feasible point x in S but not in Af*, fix) > /*; 
(ii) for all sufficiently small fik, there is an unconstrained minimizer y^ of the barrier 
function B^^{x) in strict (J-") fl int(iS), with 

B^,^{yk) = min { B^^{x) : x G strict(J') nS}. 

Thus B^^ (yk) is the smallest value of B^j,^ (x) for any x G strict(J-') fl S; 

(iii) any sequence of these unconstrained minimizers {yk} of B^^ (x) has at least one 
convergent subsequence; 

(iv) the limit point Xoo of any convergent subsequence {xk} of the unconstrained min- 
imizers {yk} defined in (ii) lies in Af*; 

(v) for the convergent subsequences {xk} of part (iv), 

lim f{xk) = f* = lim B^^{xk). 

Proof See [6]. D 

Applied to problem (3.1), the classical barrier method can be viewed as a path-following 
method that defines a path to both an optional x* and associated Lagrange multipliers A*. 
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The following theorem (stated and proved in [6] and based on results from [5, 12, 13]) 
summarizes the conditions under which a sequence of barrier minimizers converges to 
the solution of (3 . 1 ) . 

Theorem 3.3. Consider problem (3.1). Assume that the set of strictly feasible points is 
nonempty. Let x* be a local constrained minimizer of (3.1), Vf{x) = g{x), J{x) = Vc(x)'^, 
and A denote the set of indices of the active constraints at x*. Assume that the following 
sufficient optimality conditions hold at x* : 

(a) X* is a KKT point, i.e., there exists a nonempty set J^\ of Lagrange multipliers 
A satisfying 

Ma = {A : g{x*) = J{x*)^X, A > 0, and Ci{x*)Xi = \fi}; 

(b) there exists p such that Ja{x*)p > 0, where Ja{x*) denotes the active constraints 
at X*; and 

(c) there exists a; > such that p^H{x*, X)p > w\\p\\l for all A G A^a i^^^d all 
nonzero p satisfying g{x*Y p = OandJyi{x*)p > 0, where H{x*,X) = V^/(a;*) — 
Yl^i '^jV^Ci(x*) is the Hessian of the Lagrangian evaluated at x = x*. 

Assume that a logarithmic barrier method is applied in which ^k converges monotoni- 
cally to zero as k ^ oo. Then 

(i) there is at least one subsequence of unconstrained minimizers of the barrier func- 
tion B^^(x) converging to x*; 

(ii) let {x^} denote such a convergent subsequence. Then the sequence of barrier 
multipliers {Xk}, whose ith component is fik/ci{x'^) is bounded; 



(iii) limfc^oo Xix") = Xe M 



A- 



If, in addition, strict complementarity holds at x*, i.e., there is a vector X E Aix such 
that Aj > for all i E A, then 

(iv) A^ > 0; 

(v) for sufficiently large k, the Hessian matrix V'^B^^(x'') is positive definite; 
(vi) a unique, continuously differentiable vector function x(ii) of unconstrained min- 
imizers of B^(x) exists for positive jj, in a neighborhood of fi = 0; and 
(vii) lim^^o+ xijj,) = x*. 

Proof See [6]. D 

Consider the case when c{x) = x, and the constrained local minimizer x* is strictly 
positive. In this case, Theorem 3.3 reduces to the following corollary: 

Corollary 3.4. Consider problem (3.1) with c{x) = x. Suppose x* be a local constrained 
minimizer of (3.1). Further, assume that x* > 0. Let V/(x), J{x) and A be defined as 
in Theorem 3.3. Assume that the following sufficient optimality conditions hold at x*: 

(a) X* is a stationary point of f{x), i.e., g{x*) = 0; 

(b) V^/(x*) is positive definite. 

Assume that a logarithmic barrier method is applied in which fik converges monotoni- 
cally to zero as k ^ oo. Then 

(i) there is at least one subsequence of unconstrained minimizers of the barrier ffinc- 

tion B^^(x) converging to x*; 
(ii) let {x^} denote such a convergent subsequence. Then the sequence of barrier 

multipliers {A(x'^)}, whose ith component is fik/ci{x'') is bounded; 
(iii) limfc^oo A(x^) = 0; 
(iv) for sufficiently large k, the Hessian matrix 'W'^B^^{xk) is positive definite; 
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(v) a unique, continuously dijferentiable vector function x{jj) of unconstrained min- 

imizers of B^{x) exists for positive ji in a neighborhood of fi = 0; and 
(vi) lim^^o+ a;(/i) = x*. 

Proof This proof is based on the proof given in [6] for Theorem 3.3, modified for the 
case when c{x) = x and the constrained local minimizer x* is strictly positive. 

Since x* > 0, x* is in strict (J^), and thus, in the closure of strict (J-"). Assumptions 
(a) and (b) imply that x* is an isolated unconstrained minimizer of f{x). Thus, the 
conditions of Theorem 3.2 are met, implying that there is at least one subsequence of 
unconstrained minimizers of B^^ (x) converging to x*. This proves (i). 

Let {x''} denote such a convergent sequence, i.e., limfc_i.oo a^'^ = x*. Each x'' is an 
unconstrained minimizer of i?^^. (x) : 

m 
k\ \^'^„ /^k\\ (k\ „,i A (k\_ f^k 



V5^,(a:)=^(x'=)-yVQ(x')Ai(x'), where A,(x') 



Because c{x'') > (from Theorem 3.2, result (ii)), \i{x'') is strictly positive for any 
/ifc > 0. Since there are no active constraints at x* and x'' converges to x*, 

lim Ci{x^) = Ci{x*) > 0, and hence lim Xi{x^) = 0, 

k—^ca k—^oo 

for alH = 1, . . . , m, proving (ii) and (iii). 

As in the proof of (v) in Theorem 3.3, to determine the properties of V'^B^^{x^) as 
fc — 7- oo, we write the Hessian of the barrier function (3.2) as 

V^B,^{x') = V'fix') + - V A,(a;^) V'c,{x') + V ^y^Vc,{x'){Vc,{x')f. 

Since x'' -)■ x* and X{x'') -)■ as A; -)■ oo, lim^^oo V^-B^^(a;'') = V^f{x*). Thus, 
by assumption (b), for sufficiently large k, the Hessian V'^B^^{x'') is positive definite, 
proving (iv). 

To verify the existence of a unique, differentiable function a;(/i) for positive yU in a 
neighborhood of x{fik), we apply the implicit function theorem (see, for example, [11, 
p. 128] and [10, pp. 585-586]) to the n + 1 variables {x, /i). At (x'^, fi^), we know that 
the following system of nonlinear equations has a solution: 

$(x,/i) = g{x) 

The Jacobian of $ with respect to x is the barrier Hessian V^-B^(x), which was just 
shown to be positive definite at a; = x'^ and ft = fj^k- The implicit function theorem 
then implies that there is a locally unique, differentiable function x(/i) passing through 
x{nk) = x^ such that ^(x{fi), fi) = for all positive ^u in a neighborhood of /i^. 

Using continuation arguments, it is straightforward to show that the function x{fx) 
exists for all < fi < fik for all sufficiently large k, giving result (vi). 

Result (vi) is immediate from the local uniqueness of x{fi) and result (i), that Xk is a 
local unconstrained minimizer of the barrier function. D 

4. The Primal Barrier Energy Method 

Inequality constraints on u introduced in Section 2 can be enforced using a barrier 
method. Define 

J^{u) = J{u) — yU / ln(u) dx, 
Jn 
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where J{u) is defined as in (2.9). The Gateaux derivative is given by 

d 



J'^{u){v) = J'{u){v)-^^ 



/i / \u.{u + tv) 



n 
[{dVu) ■ Vf + k{u)v] dx 



t=Q 



+ / {cu — gN)v ds — fi u ^v dx. (4.1) 

Thus, the condition for stationarity of J^{u) is given by solving the following problem: 
Find u e Hlin) s.t. J'^{u){v) =0, Vi; G HliVt). (4.2) 

The second Gateaux derivative of J^(m) is given by 



J'l^{u){w,v) = I [{aVw) -Vv + k'{u)wv\ dx + cwv ds + fi u wv dx. 

(4.3) 
Thus, the Newton update for solving (4.2) is given by: 

Findu; G Hl^{n) s.t. J';{u){w,v) = -J'^{u){v), \/v G H^^iQ). (4.4) 

4.1. Discretization. We use a standard Galerkin finite element method to approximate 
the solution (4.4) in an A^-dimensional subspace Xh C X = Hq ^(Vl) fl [m_, n+]. Thus, 
we seek a solution Uh G Xh such that 

</;MM = o, ^veXh. (4.5) 

The Newton update is given by: 

Find WhEXhCX s.t. .J'^{uh){wh, Vh) = -J'^{uh){vh), ^Vh G X^,. (4.6) 
Let {(j)i}i be a basis Xh. Then, without loss of generality, let 

N N 

i i 

for some {oj} and {A}. It is sufficient to take the test functions Vh G Xh to be the basis 
functions {0j}^]^. Equation (4.6) is equivalent to solving the following matrix-vector 
equation: 

[A{uh) + fiM{uh)] W = - [G{uh) - fiH{uh)] , (4.7) 

where 



Aj{uh) 


= j;{uh){<p,Ai). 




(4.8) 


Mij{uh) 


= / {uhY'^ (t)j(t)i dx 
Jq 




(4.9) 


w. 


= A, 




(4.10) 


G^{uh) 


= ^;k)(0.), 




(4.11) 


H,{uh) 


Jn 




(4.12) 


The barrier term contributes an extra term to the system matrix, 


namely 





M = fi I u^'^(j)j(j)i dx. 



n 
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Thus, the barrier term adds a positive definite matrix to the original system matrix, and 
so, may be viewed as a regularization. In FETK, this integral is approximated using a 
high-accuracy quadrature rule, using a finite sum with fixed positive weights. 
The Newton update W defines a descent direction for 

0M = ^||GK)||^, (4.13) 

and thus, (f){uh) may be used as a merit function to enforce sufficient descent. 

At optimality, the solution Uh must lie in the strict interior of the feasible region, 
and thus, the Lagrange multipliers must be exactly zero. Because of this, there is no 
restriction that fi must be kept away from zero. Thus, /i may be steadily decreased, and 
in fact, may be set to zero-solving the original stationary problem. 

Algorithm 4.1 summarizes the primal barrier energy method: 



Algorithm 4.1: Primal Barrier Energy Method. 

Choose Mo > 0, yu > 0, r/ G (O, |), and 7 G (0, 1); 
Set k = 0; 

while not converged do 
Compute u{fi) to approximately solve (4.5): 
while not converged do 
Solve (4.7) to obtain wl; 
Use the "99% rule" to compute «max; 

Compute a G (0, Omax] such that 0(n*^ + awl) ^ 0(^h) + V^^'^ 4>{'^hY w\; 
ui ^ ui + awl; 
end do 
Uk+i <- M(/i); 
fi -ir- -ffi; 
k <- k + 1; 

end do 



In practice, each barrier function B^{x) does not have to be minimized to high pre- 
cision. Typically, each barrier function is considered sufficiently minimized when the 
norm of its gradient is either less than an fixed absolute tolerance or satisfies a relative 
tolerance based on B^{xq), where xq denotes the initial guess for the minimization (see, 
for example, [4]). 

5. Numerical Results 

The standard Newton method, a standard Newton method with safeguarding, and 
the primal barrier energy method was implemented using FETK (the Finite Element 
ToolKit; see [7] and http://www.FETK.org). These methods were used to solve the Ein- 
stein constraint equations on three single-hole domains, centered at the origin, with given 
boundary conditions on both the inner and outer boundary. Each tetrahedral mesh was 
generated by the GAMer component of FETK, which is a high-fidelity surface and vol- 
ume meshing tool based on standard simplex triangulation, subdivision, and smoothing 
algorithms (cf. [14, 15]). Details of the three meshes are given in Table 1. 

At the heart of each nonlinear solver is a linear solver (e.g., sparse direct solver or 
the conjugate-gradient (CG) method). For more ill-conditioned systems, the Newton 
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Table 1 . Meshes 





Mesh #1 


Mesh #2 


Mesh #3 


Inner radius 


50 


10 


1 


Outer radius 


100 


100 


100 


Vertices 


2089 


1436 


2820 


Simplicies 


9726 


7589 


15321 



equations may not be solved exactly; however, it is necessary that any step obtained 
from the linear solver of choice must be a descent direction. A simple backtracking line 
search is used to obtain a step that meets the sufficient decrease criteria in Algorithms 3.1 
and 4. 1 . Convergence is obtained when the norm of the nonlinear residual defining the 
PDE is less than a chosen tolerance of e = l.Oe-07, i.e., 

\\G{uk)\\2<e. (5.1) 

For each example, the energy barrier method initialized fi and then decreased /i when- 
ever the iterate m^ satisfied 

||/(4)||2<max{e,||/K)||2,e,} (5.2) 

where /(m^) = G{u^) — fiH{u'^), m° denotes the initial iterate after decreasing /i, and 
e^ = max{min{0.1, /i}, e}. This choice of e^ allows each subproblem to be solved to 
greater accuracy as /i is decreased. 

Reasonable choices for parameters for the Lichnerov equation with boundary condi- 
tions include those given in Examples 1-2, in Section 2.2. The first two examples are 
with these choices of parameters. 

Example 1. Set a = 1.0, R = 1.0, r = 0.1, a = 0.2, p = 0.1, c = 1.0, and qn = -1-0. 
We pick n to have Robin boundary conditions on both the inner and outer boundary. 
The presence of the negative exponents acts as a natural barrier function, preventing in- 
feasible iterates. For this reason, to obtain a positive solution it is sufficient to add a 
safeguarding procedure such as the so-called "99% rule" (see Section 3) to the standard 
Newton method. 

Table 2 gives the results of using the standard Newton method, Newton's method with 
safeguarding using the "99% rule", and the primal energy method to solve the Lich- 
nerovich equation. For each solver, the initial guess was a vector of all ones. We list 
the number of iterations ("itns"), residual, and the signs of the entries in the vector of 
coefficients for Uh (i.e., + denotes all the entries are strictly positive, — denotes all the 
entries are strictly negative, and +/— denotes both positive and negative entries). Note 
that there is one linear solve per iteration, and thus, the number of iterations is also the 
number of linear solves required by each method. All four solvers converged on all three 
meshes to strictly positive solutions. 

In Table 2, the energy barrier method is reported with two different initial values of /i. 
First, the energy barrier method was run with fi = 0.0, making it numerically equivalent 
to Newton's method with safeguarding. (For nonconvex problems, we do not expect 
the energy barrier method to converge to a strictly positive solution with this choice of 
yu). For illustrative purposes, the results with /i = 1.0 are also reported; for this test, 
when each subproblem was sufficiently solved (i.e., (5.2) was satisfied), we reduced /i 
by a factor of 1/10. However, for this convex problem, allowing faster reductions of 
H leads to fewer overall iterations. In fact, reducing ^u by a factor of 1/100 led to 14 
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iterations on mesh #1, 16 iterations on mesh #2, and 16 iterations on mesh #3. For more 
difficult problems, reducing jj, too quickly will impede convergence. Even though faster 
reductions in fi would lead to results more similar to the Newton methods, for the results 
in Table 2, we chose to display results with a reasonable, commonly accepted reduction 
of 1/10. Also, it is worth noting that the different residual values for the algorithms are 
inconsequential in that convergence only requires (5.1) to be satisfied. (In these cases, 
Newton's method was fortunate in that quadratic convergence led to a much smaller 
residual for each final iterate.) 

In Table 2 we see that a positive solution was obtained by all methods using an initial 
guess of all ones. As previously noted, a barrier method is not required to obtain a 
strictly positive solution since the energy functional has a natural barrier in the form of 
negative coefficients of u. Thus, we expect that a safeguarded Newton method will be 
sufficient to recover a positive solution. In this example, Newton's method with and 
without safeguarding produced the same iterates. It is of interest to note that a strictly 
negative solution can be recovered on all three meshes by using the standard Newton 
method together with the initial guess of a vector of all negative ones. 

Table 2. Example #i. 





Mesh #1 


Mesh #2 


Mesh #3 


method 


itns 


resid 


sign 


itns 


resid 


sign 


itns 


resid 


sign 


Newton (standard) 


6 


2.71e-12 


+ 


6 


5.15e-12 


+ 


6 


4.73e-12 


+ 


Newton (safeguai-ded) 


6 


2.71e-12 


+ 


6 


5.15e-12 


+ 


6 


4.73e-12 


+ 


Barrier energy (//q — 0.0) 


6 


2.71e-12 


+ 


6 


5.15e-12 


+ 


6 


4.73e-12 


+ 


Barrier energy (//q = 10) 


22 


6.15e-08 


+ 


24 


1.19e-08 


+ 


24 


1.17e-08 


+ 



Example 2. Let a = 2, i? = -1000, r = V72, a = V48, and p = l/vr. For this 
example, the inner and outer boundaries of all three meshes have a Robin boundary 
condition with c = 2 and g^ = 10. 

On the first mesh, the standard Newton method failed to converge in 50 iterations, 
denoted by the asterisk in Table 3; in fact, at the 100th iteration, the current approxi- 
mation to Uh contained both positive and negative entries-the standard Newton method 
was unable to maintain a strictly positive solution. When the initial guess was set to 
a vector whose entries were all 10, Newton's method converged to a positive solution; 
when the initial guess was set to a vector whose entries were all -10, the standard New- 
ton method returned a strictly negative solution; and when the initial guess was set to a 
vector whose enries were all -1, the standard Newton method returned a strictly positive 
solution. (Newton's method without safeguarding was unpredictable-one could recover 
a positive solution even when starting with an initial negative guess.) Meanwhile, with an 
initial guess of all ones, the safeguarded Newton method converged to a strictly positive 
solution. The barrier energy method converged with fiQ = I and subsequent reductions 
in /i of 1/10, as in Example 1. 

On the second mesh, all methods converged to the same strictly positive solution when 
the initial guess was all ones. (The barrier energy method was run with jiq = 1 and 
subsequent reductions in /i of 1/10, as in Example 1.) A strictly negative solution can 
be recovered by using the standard Newton method together with the initial guess of a 
vector whose entries are all -10. 



BARRIER METHODS FOR CRITICAL EXPONENT PROBLEMS 



15 



On the third mesh, Newton's method without safeguarding did not converge within 
the first 100 iterations and had both positive and negative components. A negative so- 
lution was obtained by the standard Newton method by starting with a vector whose 
entries were all -10. Meanwhile, Newton's method with safeguarding made insignificant 
progress after the seventh iteration, and thus, failed to converge in 100 iterations. The 
safeguarded method failed because the initial safeguarded Newton steps took some of 
the coefficients of u^ very close to the boundary and the following Newton directions 
continued to point in the direction of negative numbers for these components. In this 
event, the "99% rule" allows subsequent steps of only negligible size along the Newton 
direction at each iteration. As a result, the safeguarded Newton method is unable to make 
any real progress each subsequent iteration, and thus, fails to converge. 

In this example, we see that even though there is a natural barrier in the Lichnerovich 
equation, there are additional benefits that a barrier function approach can offer: By 
setting the parameter ^ to be large enough, we can alter the pure Newton direction, pre- 
venting the initial Newton iterates from getting too close to the boundary. For example, 
setting fiQ to be 10, 20, 30, or 40, the energy method's iterates took large initial steps to 
the boundary, preventing convergence; however, with fiQ = 50, the algorithm converged 
to a positive solution. (The values in Table 3 are with /xq = 50 and subsequent reductions 
in /i of 1/10.) 

The energy barrier method was the only method to find a strictly positive solution to 
this problem on all three meshes. 

Table 3. Example #2. 





Mesh #1 


Mesh #2 


Mesh #3 


method 


itns 


resid 


sign 


itns 


resid 


sign 


itns 


resid 


sign 


Newton (standard) 


* 


4.69e+19 


+/- 


11 


9.80e-08 


+/- 


* 


4.92e+20 


+/- 


Newton (safeguarded) 


9 


2.89e-09 


+ 


7 


1.37e-08 


+ 


* 


2.01e+07 


+ 


Barrier energy 


16 


4.02e-08 


+ 


16 


7.67e-08 


+ 


17 


9.47e-09 


+ 



Example 3. Consider following Yamabe problem: 



-8Am 



0, 



where p(r) = 1/r'^, and r is the Euclidean distance. (This choice of p(r) was motivated 
by equation (41) in [1]). For this example, we impose the Dirichlet condition u = 1 
on both the inner and outer boundaries of all three meshes. For this example, the initial 
guess was taken to be a vector of ones; given the Dirichlet condition, this is a reasonable 
starting point. 

Table 4 reports the results of each solver on this problem. Neither backtracking nor a 
barrier method was required to solve this problem. 



Table 4. Example #3. 





Mesh #1 


Mesh #2 


Mesh #3 


method 


itns 


resid 


sign 


itns 


resid 


sign 


itns 


resid 


sign 


Newton (standard) 


1 


1.91e-08 


+ 


2 


1.16e-12 


+ 


3 


1.18e-12 


+ 


Newton (safeguarded) 


1 


1.91e-08 


+ 


2 


1.16e-12 


+ 


3 


1.18e-12 


+ 


Barrier energy (^=1.0) 


18 


1.13e-12 


+ 


22 


1.25e-12 


+ 


23 


1.34e-12 


+ 
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Example 4. For this example, we modify the problem in Example 3 to include the extra 
term in (1.1): 

—8Au u = —p{r)u^, 

8 

where p(r) = 1/r^, and r is the Euclidean distance. Also, assume Dirichlet boundary 
conditions of n = 1 on both the inner and outer boundaries. For this example, the initial 
guess was taken to be a vector of ones. 

Table 5 contains the results on all three meshes. The standard Newton method con- 
verged quickly to a solution with positive and negative components on each mesh. New- 
ton's method with safeguarding did not converge on any mesh. On the all three meshes, 
the method took initial large steps to the boundary and made negligible progress after ten 
iterations. The primal barrier energy method with fiQ = 1 converged on the first mesh, 
but this initial value of /i on the second mesh was too small. With fiQ = 10, the primal 
barrier energy method converged to a strictly positive solution on all three meshes. 

Tables. Example #4. 





Mesh #1 


Mesh #2 


Mesh #3 


method 


itns 


resid 


sign 


itns 


resid 


sign 


itns 


resid 


sign 


Newton (standard) 


11 


2.71e-08 


+/- 


16 


9.65e-09 


+/- 


18 


8.97e-ll 


+/- 


Newton (safeguarded) 


* 


9.16e+03 


+ 


* 


1.64e+04 


+ 


* 


1.47e+04 


+ 


Barrier energy (/i=10.0) 


17 


1.94e-ll 


+ 


18 


2.86e-ll 


+ 


18 


2.69e-ll 


+ 



6. Conclusion 

In this article we considered both the design and the analysis of a certain class of non- 
convex optimization-based numerical methods for approximating positive solutions to 
nonlinear geometric elliptic partial differential equations containing critical exponents. 
As noted, these types of problems arise regularly in geometric analysis and mathemati- 
cal physics; our primary interest here was Yamabe problem and the Einstein constraint 
equations. The difficulty one faces with these problems are the simultaneous presence of 
several challenging features, including spatial dimension n ^ 3, varying and potentially 
non-smooth coefficients, critical (or even super-critical) nonlinearity, non-monotone non- 
linearity (arising from a non-convex energy), and spatial domains that are typically Rie- 
mannian manifolds rather than simply open sets in i?". For these types of problems, there 
may be multiple solutions, although only positive solutions typically have mathemati- 
cal and physical meaning. This creates additional complexities in both the theory and 
numerical treatment of such problems, as this feature introduces both non-uniqueness 
as well as the need to incorporate an inequality constraint into the formulation. As 
a practical approach for treating these difficulties, we considered numerical methods 
based on Galerkin-type discretization, covering any standard bases construction (finite 
element, spectral, or wavelet), and the combination of a barrier method for nonconvex 
optimization and global inexact Newton-type methods for dealing with nonconvexity 
and the presence of inequality constraints. After giving an overview of barrier methods 
in non-convex optimization, we then developed and analyzed a primal barrier energy 
method. We then presented a sequence of numerical experiments using this type of bar- 
rier method, based on a particular Galerkin method, namely the piecewise linear finite 
element method, leverage the FETK modeling package. In the experiments, the negative 
pole in the Hamilitonian constraint provided a "natural" barrier, aiding the convergence 
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of Newton methods. In this setting, a barrier method will often be unnecessary; how- 
ever, in some cases the numerical experiments showed that allowing for a flexible barrier 
parameter can be helpful (see Example 2, Mesh #3). The experiments also confirmed 
that on some classes of Yamabe problems, a solution could not be found without the use 
of a barrier method (see Example #4), suggesting that barrier methods are more useful 
on critical exponent problems without singularities that arise from lower-order negative 
exponent terms. 

Although we considered here only scalar elliptic equations with variational structure, 
that is, they that arise as the Euler condition for stationarity of an underlying (usually 
nonconvex) energy, more generally these types of critical exponent problems may arise as 
part of a more complex elliptic system. A prime example is the coupled Hamiltonian and 
momentum constraints in the Einstein equations [8, 9]. While the Hamiltonian constraint 
(as well as the momentum constraint) alone has variational structure, when combined as 
a system there is in fact no variational structure to exploit. Nevertheless, the ideas in this 
paper can be applied by using alternative formulations of the Einstein constraints, and 
will be pursued in a second article. 
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